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Abstract 

We  develop  a  formulation  for  the  single-fluid / two-temperature  equa¬ 
tions  for  simulating  two-species,  compressible,  non-equilibrium  plasma 
flows.  The  divergence-free  condition  of  the  magnetic  field  is  enforced 
via  the  characteristic  decomposition  of  an  extended  nine- wave  system. 
The  source  terms  are  modified  appropriately  to  improve  energy  and 
momentum  conservation  accuracy.  A  spectral//ip  element  algorithm  is 
employed  in  the  discretization  combined  with  a  discontinuous  Galerkin 
formulation  for  the  advective  and  diffusive  contributions.  The  for¬ 
mulation  is  conservative,  and  nronotonicity  is  enforced  by  appropri¬ 
ately  lowering  the  spectral  order  around  discontinuities.  A  new  MHD 
flux  introduced  here  is  the  MHD-HLLC  (Harten-Lax-van  Leer  Contact 
wave)  flux  that  preserves  nronotonicity  and  resolves  contact  disconti¬ 
nuity  better.  Exponential  convergence  is  demonstrated  for  a  magneto- 
hydrostatic  problem.  Two  tests  are  presented  using  the  new  MHD- 
HLLC  flux.  Also,  the  differences  between  the  single-temperature  and 
the  two-temperature  models  are  presented  for  two-dimensional  plasma 
flows  around  bluff  bodies  are  simulated. 
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1  Introduction 


Plasmas  can  be  modelled  accurately  using  kinetic  theory,  especially  partially 
ionized  plasmas.  However,  this  involves  solutions  of  the  seven-dimensional 
Boltzmann  equation  coupled  with  Maxwell’s  equations,  which  is  prohibitively 
expensive.  Particle  based  methods,  such  as  DSMC,  are  possible  alternatives 
but  for  efficiency  they  need  to  be  coupled  with  continuum  fluid  equations. 
Such  hybrid  methodologies  have  been  used  successfully  recently  in  the  sim¬ 
ulation  of  ion  thruster  plumes,  but  hybrid  kinetic-continuum  methods  are 
still  under  development  [1,  2],  and  open  issues  with  the  DSMC  remain  the 
treatment  of  electrons  as  well  as  the  modelling  of  charged  particle  collisions. 

Continuum-based,  i.e.  purely  fluid  approaches,  have  been  successful  in 
describing  the  macroscopic  features  of  high  density  plasmas  in  many  diverse 
applications  [3,  4,  5,  6,  7,  8,  9,  10].  They  are  derived  from  the  Boltzmann 
equation  by  taking  appropriate  moments  for  each  species.  The  standard 
mathematical  description  is  that  of  single-fluid  MHD  with  magnetic  and  gas 
dynamic  viscous  effects.  However,  a  single-fluid  MHD  description  has  its  lim¬ 
itations  as  it  cannot  account  for  local  thermodynamic  non-equilibrium  effects 
and  cannot  consider  non-neutral  regions  and  sheath  interactions.  To  this  end, 
two- fluid  plasma  models  [11]  and  corresponding  solvers  have  been  under  de¬ 
velopment  more  recently  [12].  They  can  overcome  certain  limitations  of  the 
single-fluid  MHD  model  such  as  the  Hall  effect  and  diamagnetic  terms,  which 
model  contributions  to  ion  current  and  the  finite  Larmor  radius  of  the  plasma 
constituents.  However,  they  still  assume  local  thermodynamic  equilibrium 
within  each  fluid.  From  the  computational  standpoint,  the  two- fluid  model 
is  much  more  complex  system  to  solve,  especially  for  large  values  of  the  Hall 
parameter,  and  approximate  Riemann  solvers  are  still  under  development 

I12]. 

In  between  the  single- fluid  and  the  two- fluid  models  for  plasmas  is  the  two- 
temperature  model.  It  can  account  partially  for  the  energy  transfer  between 
heavy  species  and  electrons,  and  it  is  computationally  more  tractable.  There 
can  be  many  applications  for  which  the  electron  temperature  differs  from  the 
heavy  particle  temperature  [13,  14,  15].  For  example,  experiments  in  [13] 
with  boundary  layers  of  pure  NaK  seeded  argon  showed  that  the  electron 
temperature  was  considerably  higher  than  the  gas  temperature.  Moreover, 
the  electrical  conductivity  and  other  transport  coefficients  in  the  conservation 
laws  depend  strongly  on  the  electron  temperature.  This  is  also  true  for 
emerging  applications  in  microfluidics,  such  as  micro-pulsed  plasma  thrusters 
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(micro-PPTs)  where  viscous  effects  dominate,  and  the  flow  may  be  in  a 
transitional  state. 

The  objective  of  the  current  work  is  to  develop  robust  high-order  al¬ 
gorithms  for  a  single-fluid/two-temperature  plasma  extending  our  previous 
work  on  the  standard  MHD  model  [16].  The  use  of  high-order  accuracy 
addresses  effectively  the  small-scale  requirements  of  compressible  MHD  tur¬ 
bulence  [17],  as  well  as  the  extra  resolution  required  in  long-time  integration 
[18].  In  addition  to  our  work,  other  recent  efforts  to  develop  effective  high- 
order  methods  for  plasma  flows  have  been  reported  in  [19,  20,  21]. 

Discontinuous  Galerkin  methods  [22]  address  two  of  the  main  difficulties 
in  employing  high-order  discretization  for  the  solution  of  hyperbolic  conser¬ 
vation  laws  are: 

1.  Maintaining  monotonicity  for  non-smooth  solutions,  and 

2.  preserving  conservativity. 

In  the  MHD  framework,  such  difficulties  are  compounded  by  the  impo¬ 
sition  of  the  divergence- free  condition  for  the  magnetic  held,  which  results 
in  a  loss  of  the  hyperbolicity  of  the  ideal  MHD  equations.  This  condition 
has  been  dealt  with  by  employing  staggered  grids  in  the  work  of  Evans  & 
Hawley  [23],  which  was  extended  by  Peterkin  et  al.  [5].  However,  such  an 
approach  cannot  be  easily  incorporated  in  high-order  discretizations.  Al¬ 
ternative  approaches  include  the  operator-splitting  algorithm  proposed  by 
Zachary,  Malagoli  &  Colclla  [4]  and  the  the  development  of  extended  Rie- 
mann  solvers  by  Powell  [6];  the  latter  is  easily  extended  to  multi-dimensions 
and  also  to  high-order  discretization.  In  some  approaches,  the  divergence- 
free  condition  is  not  imposed  directly  during  time-stepping  but  the  initial 
conditions  are  projected  in  the  divergence-free  space.  Assuming  that  the  flux 
of  divergence  of  the  magnetic  held  satisfies  a  homogeneous  discrete  parabolic 
equation  with  homogeneous  boundary  conditions,  then  this  will  lead  to  zero 
discrete  divergence  at  all  times.  However,  in  practice,  discretization  errors 
or  other  inconsistencies  may  trigger  large  divergence  errors  for  such  cases. 

In  this  paper  we  extend  the  approximate  Riemman  solver  of  Powell  [6] 
to  a  nine-wave  system  to  account  for  the  divergence-free  condition  and  the 
extra  electron  energy  equation.  Compared  to  the  Riemann  solver  [6]  and 
also  our  previous  work  in  [16],  here  we  modify  the  source  terms  to  preserve 
conservativity  and  enhance  accuracy.  We  employ  a  spectral/hp  element  dis¬ 
cretization  [24]  based  on  tensor-product  polymorphic  elements  for  the  mesh 
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macro-skeleton.  The  spectral  order  of  the  Jacobi  polynomials  in  the  trial  ba¬ 
sis  is  variable  in  order  to  accommodate  different  solution  requirements,  i.e. 
low-order  for  shock  capturing  or  high-order  boundary  layer  resolution. 

To  preserve  pressure  positivity,  a  MHD-HLLC  flux  [25]  is  implemented  in 
the  discontinuous  Galerkin  method  for  solving  compressible  plasma  flows.  In 
addition  to  capturing  the  effects  of  contact  waves,  the  MHD-HLLC  flux  also 
resolves  Alfven  and  slow  waves  better  than  the  HLL  (Harten-Lax-van  Leer) 
flux  and  the  Lax-Friedrichs  flux. 

As  a  model  problem  we  use  plasma  flow  past  a  cylinder.  Depending 
on  the  specific  conditions,  we  find  that  the  electron  temperature  can  be 
substantially  different  than  the  temperature  of  ions,  and  correspondingly 
this  may  affect  the  velocity  field.  To  appreciate  the  differences  we  compare 
the  two-temperature  model  with  the  standard  MHD  model  under  the  same 
wall  thermal  condition. 

The  paper  is  organized  as  follows:  In  section  2  we  present  the  formulation 
and  briefly  summarize  details  of  the  implementation.  In  section  3,  we  first 
test  the  convergence  rate  of  the  algorithm  for  an  analytical  problem.  We  then 
present  tests  using  the  new  MHD-HLLC  flux,  and  subsequently  we  simulate 
plasma  flow  past  a  cylinder  in  the  subsonic  and  supersonic  regimes.  Finally, 
we  conclude  in  section  4  with  a  few  remarks. 

2  TWO-TEMPERATURE  PLASMA  EQUATIONS 
2.1  Governing  Equations 

The  non-dimensional  governing  equations  for  single-fluid/two-temperature 
plasma  for  compressible  magneto- hydrodynamics  (MHD)  can  be  expressed 
in  conservative  form  as  (see  derivation  in  Appendix  B): 

1.  Mass  Conservation : 


dp 

dt 


-V  .  (pv) 


2.  Momentum  Conservation : 


(1) 


d(pv) 

dt 


— V  •  (pvv* 


BB*  +  (p  +  -|B|2)/  —  — t, 

Z  Oni 


(2) 
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3.  Magnetic  Field : 


(3) 


<9B  1 

-  =  -Vx(  Bxvt-VxB| 


4.  Total  Energy  Conservation : 

~  =  —  V  ■  [{Etot  +  p)v  +  (^|B|2I  —  BB')  •  v  —  v  •  Tj  —  77— v  •  re 

eye  z  O915 


-J-(B^VB-V(i|B|2))--^  ... 

Ej f'  ^  O'ljp-L  /  c 


-  VTP 
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c  ,pr. 

Uyi-1  I  % 


VTA 


(4) 


5.  Electron  Internal  Energy  Conservation : 

t  -  -v'i'£«+^-sikvrJ 

-re  :  Vv  +  — (V  x  B)  •  (V  x  B) 


v  •  Vpe 


6.  Magnetic  Flux  Constraint: 


Sv 


Sr 


(5) 


V  •  B  =  0 


(6) 


7.  Ohm's  Law : 


E  =  ijJ-vxB 


(7) 


Here  we  define: 


Etot  =  7  P  .  +  x(Pv  •  v  +  B  •  B),  p  =  pi+pe ,  ee  =  — (8) 
(7  -  1)  2  7-1 


The  stress  tensor  for  ions  and  electrons  is  defined  as: 

2  2 
Ti  =  (djVH  +  S(Vy)  -  - V  •  v/i.j.  r,  =  (djVei  +  v.j  -  -V  ■  v./i,, 


(9) 


All  other  parameters  are  as  defined  in  table  1.  The  subscript  ‘i’  denotes 
ions  while  the  subscript  ‘e’  denotes  electrons. 
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The  above  Ohm’s  law  is  simplified  from  the  generalized  Ohm’s  law  in  the 
limit  of  small  Larmor  radius  approximation.  The  generalized  Ohm’s  law  can 
be  expressed  as: 


E 


m 


+  v  x  B  + 


j  x  B 

ne 


Vpe 

ne 


(10) 


and  thus 


j  x  B/ne  Vpe/ne  rLi 
-  - 

v  x  B  v  x  B  L 


(ii) 


where  rjA  is  the  ion  Larmor  radius  and  L  is  the  scale-length  of  the  fluid 
motion.  We  conclude  that  the  second  and  third  terms  on  the  right-hand- 
side  of  equation  (10)  can  be  neglected  if  the  ion  Larmor  radius  is  very  small 
compared  to  the  characteristic  length  scale  of  the  fluid  motion,  i.e.,  ru/ L  « 
1.  Specifically,  when  we  consider  the  length  scale  of  the  fluid  motion  to  be 
very  small  and  close  to  the  ion  Larmor  radius,  then  we  have  to  include  the 
two  additional  terms  in  equation  (10)  and  use  the  generalized  Ohm’s  law 
instead. 


Alternatively,  in  flux  form,  the  above  conservation  equations  can  be  expressed 
compactly  as 


d»U 

~dt 

+ 

where  all  flux  and  source  terms  are  defined  in  detail  in  the  Appendix  A.  The 
state  vector  is  defined  as:  U  =  (p,  pu,  pv,  pw,  Bx,  By,  Bz,  Etot,  ee) 

2.2  The  V  •  B  =  0  Constraint 

The  presence  of  the  V  •  B  =  0  constraint  implies  that  the  equations  do  not 
have  a  strictly  hyperbolic  character.  It  has  been  shown  in  [26]  that  even  a 
small  divergence  in  the  magnetic  fields  can  dramatically  change  the  character 
of  results  from  numerical  simulations.  In  our  work,  we  adopt  an  approach 
which  was  developed  originally  by  Powell  in  [6].  The  idea  is  to  re-formulate 
the  Jacobian  matrix  to  include  a  “ninth- wave” ,  i.e.,  the  divergent  mode  that 
corresponds  to  velocity  u.  This  way  the  degeneracy  associated  with  the 
divergence-free  condition  is  avoided  while  the  rest  of  the  eigenvalues  of  the 
Jacobian  remain  the  same. 


OF 


Ideal 


Qjpldeal 


dF 


Ideal 


dx 

dFlisc 

dx 


dz 


dy 

dFvvisc  dFYisc  „ 

H - ^ - 1 - ^ - 1-  Smhd 


dy 


dz 


(12) 

(13) 
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The  primitive  Jacobian  matrix  Ap  for  single- fluid/ one-temperature  equa¬ 
tions  has  the  form,  in  three-dimensions, 

0  0  ' 

B,  i 

p  p 

0  1 
p 

—  0 

p 

0  0 
0  0 
u  0 
0  u 

Considering  that  p  =  Pi  +  pe,  then  the  primitive  Jacobian  matrix  Ap  for 
single- fluid/two-temperature  equations  in  three-dimensions  can  be  described, 

i.e. 


Ap  — 


u 

0 

0 

0 

0 


p 

u 

0 

0 

0 


0  By 

0  Bz 

0  7 p 


0 

0 

u 

0 

0 

~BX 

0 

0 


0 

0 

0 

u 

0 

0 

-  B , 


0 

Bx 


0 


—v 

—w 

-(7  —  l)u  •  B 


0 

Bjt 

p 

Bx 

P 

0 

0 

u 

0 

0 


Ap  — 


up  0  0 

0  u  0  0 

0  0  u  0 

0  0  0  u 

0  0  0  0 

0  By  -Bx  0 

0  Bz  0  -Bx 

0  7p  0  0 

0  7Pe  0  0 


0 

_ Bx_ 


P 

0 


—v 

—w 

—  (7  —  l)u  •  B 
0 


0  0  0  0  ■ 

By  Bz  1  Q 

p  p  p 

^0^0 
p  n  p 

0^00 

0  0  0  0 

u  0  0  0 

0  u  0  0 

0  0  u  0 

0  0  0  m 


To  modify  the  governing  equations  so  as  to  make  Ap  non-singular,  using 
Powell’s  criteria  presented  in  [6],  Ap  is  modified  to  be  Ap: 


AL 


up  0 
0  u  0 

0  0  u 

0  0  0 

0  0  0 

0  By  —Bx 

0  Bz  0 

0  7p  0 

0  7  pe  0 


0  0  0 

0  0^ 
0  0^ 
u  0  0 

0  u  0 

0  0m 

-Bx  0  0 

0  0  0 

0  0  0 


0  0  0  1 


„  p 

^00 

0  0  0 

0  0  0 

M  0  0 

0  M  0 

0  0m 
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This  modification  effectively  corresponds  to  adding  a  source  term  propor¬ 
tional  to  V  •  B, 

& Poweii  =  -(V  •  B)(0,  Bx,  By,  Bz,u,v,w,u  •  B,  0)T 

to  the  right-hand-side  of  all  evolution  equations.  We  note  that  this  source 
term  does  not  shift  the  physical  solution  since  V  •  B  is  imposed.  However, 
some  local  accumulations  may  occur,  especially  at  stagnation  points  for  which 
v  =  0.  In  these  cases,  it  may  be  necessary  to  add  the  HclmHoltz  projection, 
see  [26]. 

Next,  we  show  how  we  need  to  modify  these  source  terms  in  order  to 
better  maintain  pressure  positivity.  In  solving  the  MHD  system,  the  pressure 
is  a  derived  variable.  Specifically,  it  is  obtained  by  subtracting  off  the  kinetic 
energy  and  magnetic  energy  from  the  total  energy.  However,  in  applications 
of  micropropulsion,  magneto-spheric  physics  and  astrophysics,  the  pressure 
can  be  several  orders  of  magnitude  smaller  than  either  the  kinetic  energy 
or  the  magnetic  energy.  Thus,  small  discretization  errors  in  the  total  energy 
can  produce  situations  where  the  pressure  might  become  negative.  This  leads 
to  an  unacceptable  physical  situation.  As  long  as  the  regions  in  front  of  a 
magneto-sonic  shock  have  positive  pressure,  negative  pressures  would  not  be 
produced  in  magneto-sonic  shocks. 

Janhunen  [27]  has  reported  that  the  solution  of  the  Riemann  problem  for 
Powell’s  equations  for  left-  and  right-states  with  positive  fluid  pressures  may 
contain  unphysical  intermediate  state  with  negative  fluid  pressure,  pressure 
positivity,  as  well  as  energy  and  momentum  conservation  could  be  regained 
by  discarding  the  source  terms  in  the  energy  and  momentum  equations,  so 
that  the  source  term  proportional  to  V  ■  B  becomes: 

S  =  —(V  •  B)(0,  0,  0,  0,  u,  v,  w,  0,  Of, 

2.3  Implementation  of  the  Inviscid  Terms 

We  evaluate  the  inviscid  fluxes  and  their  derivatives  in  the  interior  of  the 
elements  and  add  correction  terms  (jumps)  for  the  discontinuities  in  the 
flux  between  any  two  adjacent  elements.  In  order  to  evaluate  the  Euler 
flux  at  an  element  interface,  we  use  an  one-dimensional  Riemann  solver  to 
supply  an  upwinded  flux  there,  see  below.  At  a  domain  boundary,  we  provide 
far  field  conditions  and  treat  the  exterior  boundary  as  the  boundary  of  a 


“ghost”  element.  This  way  we  can  use  the  same  Riemann  solver  at  all  element 
boundaries. 

We  linearize  the  one  dimensional  flux  ~pIdeal  in  the  normal  direction  to 
a  shared  element  boundary  using  the  average  of  the  state  vector  at  either 
side  of  the  element  boundary.  That  is,  since  F Ixdeal  is  a  nonlinear  function  of 
the  state  vector,  we  use  the  average  state  to  form  an  approximation  to  the 
Jacobian  of  the  flux  vector  Ac. 

The  Jacobian  matrix  for  the  flux  vector  of  the  evolution  equations  ex¬ 
pressed  in  primitive  variables  is  simpler  than  in  the  conserved  form.  Thus, 
we  will  perform  the  linearization  about  the  primitive  form  and  transform  to 
the  conserved  form.  The  left  and  right  eigenvectors  of  the  primitive  Jacobian 
matrix  Ap,  similar  to  the  results  shown  in  [6],  are: 

Entropy  wave : 

Ae  u 

le  =  (1,0,  0,0,  0,0,0, -4,0) 

gt 

'r  e  =  (1,0,0,0,0,0,0,0,0)4 

Alfven  Waves : 


Bx 

A  a  =  U±  — — 

Vp 

L  =  ^{0A~Pz,PyA±^,T^zA0) 
ra  =  ^(0,  0,  py,  0,  ±&VP,  Tf3yVP,  0,  0)T 

Fast  waves: 


A/  =  u±cf 

l  1  /n  .  Q  a  Q  Q  n 

If  =  X  9(0,  i.OCfCf,  T OisCsPxjjy ,  -f-OlsCsfjxPz 1  0,  —  ,  —  ,  ,  0) 

2a2  ^TP  \fP  P 

77  =  ( paf ,  ±afcf,  TotsCsfixPy,  Tascs(3x(3z,  0,  asf3ya^/p,  asf3za^/p,  afjp,  appe 


Slow  waves: 


Xs  =  uics 
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ls  =  ^^(0,^oiscs,±afcfPxPy,±afcfPxpz,0,:^-  ^  ,—  ^Z,~, 0 ) 

=  (pa„,  ±aaca,  ±afCf/3xPy,  ±afcfpxpz,  0,  - afpya y/p,  -afpza^/p,  as^p,  as^pPf 

Compared  to  [6],  we  have  an  extra  wave  which  corresponds  to  the  electron 
energy: 


lee  =  (0,0,  0,0,  0,0,  0,0,1) 


r,e  = 


IPe 


,0,0,  0,0,  0,0, 0,1)* 


Here: 


^  *  \  2 


7P  +  B  •  B 

P 

1 


cl  =  -  I  (a*)2  +  \l (a*)4  -  4  —  ~  |  ,  c;  =  -  |  (a 


,7 PBx\ 


*  \  2 


p2  r  s  2 


,2  „2 


r>2  «2 


af  - 


a"  —  ct  o  Ci  n 

- -  =  — - 

r2  „2  ’  s  «2  _2 

C/  C*  C/  Cs 


Px  =  sgn(Bx),  Py  = 


B„ 


Bl  +  52 


,  Pz  = 


B, 


Bl  +  52 


7*^4 


,7ggTl 

P2  ) 


We  can  transform  between  the  primitive  variables  W  and  conserved  variables 
U  with  the  following  transform: 

A  _  dU  A  <9W 

c  “  aw-  p ~du 

where 

U  =  (p,  pu,  pv,  pw,  Bx,  By ,  Bz,  Etot,  pee) 
are  the  conserved  variables ,  and 


W  =  (p,u,v,w,  Bx,  By,  Bz,p,pe) 
are  the  primitive  variables.  This  gives: 
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dU 

dW 


1  0  0 

up  0 
v  0  p 

w  0  0 

0  0  0 

0  0  0 

0  0  0 

^  pu  pv 

0  0  0 


0  0  0 

0  0  0 

0  0  0 

p  0  0 

0  1  0 

0  0  1 

0  0  0 

pw  Bx  By 

0  0  0 


and 


dW 

~dV 


10  0  0 

i  0  0 

p  p 

-^0^0 

p  p 

0  0  i 

p  p 

0  0  0  0 

0  0  0  0 

0  0  0  0 


0 

0 

0 

0 

1 

0 

0 


•  u  — ju  —'jv  —7 w  —7 Bx 

0  0  0  0  0 


0  0  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 

1  0  0 

Bz  — l—  0 

0  0^ 

7-1 


0  0 

0  0 

0  0 

0  0 

0  0 

1  0 

0  1 

-7  By  -7  Bz 

0  0 


0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
0  0 
7  0 
0  7 


where  7  =  7  —  !. 


Figure  1:  Interface  conditions  between  two  adjacent  triangles. 

We  are  now  in  a  position  to  evaluate  different  fluxes  at  the  element  bound¬ 
aries.  The  formulations  of  different  fluxes  we  employ  are: 
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1.  Flux  employed  in  [28]: 

1  dll  k= 9 

F^IU^Ue)  =  -[F(U7)  +  F(Ub)-^  £afc|Afc|rk]  (14) 

dW 

ak  =  Ik’ -^(Ub-U/), 

2.  Lax-Friedrichs  flux: 

F  Lax— Friedrichs  (U7,UE)  =  ^[F(Uj)  +  F(Ue)  —  Dmax(UE  —  U/)]  (15) 

where  Dmax  =  max(\Xk\)  -  maximal  absolute  value  of  the  eigenvalues, 
‘I’  denotes  interior  and  ‘E’  denotes  exterior  of  the  element  (see  figure 
1).  Here,  the  lfc  and  rk  are  the  ordered  left  and  right  eigenvectors  of 
the  primitive  Jacobian  matrix.  We  have  to  apply  the  operator  to 
the  right  eigenvectors  to  calculate  the  conserved  flux.  The  \'ks  are  the 
wave  speeds  associated  with  the  eigenvectors. 

A  new  MHD  flux  introduced  here  is  based  on  the  MHD-HLLC  flux 
presented  in  [25]  that  can  preserve  positivity  and  improve  resolution, 
especially  at  contact  interfaces. 

3.  MHD-HLLC  interface  flux: 


F_ffLLc(U/,  Us)  —  < 


F,  if  Sr>0 


F*  if  5/  <  0  <  S M 


F*e  if  SM  <  0  <  SE 


Fr  if  5r  <  0 


where  Fj  and  F^  are  defined  as: 


p*  =  p*  _  Si(SE  -  gM)AU, 


FL  =  F*  + 


(SE  ~  Sj) 

SE(SM  -  Sj) 
(SE  -  Sj) 


AU* 


(16) 


(17) 
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with  F*  =  SeFi-SjFe+SiSeCUe-U!)  and  AU*  =  u*,  _  XJ}.  The  wave-speeds 
Si,  Sm  and  Se  are  defined  as: 

'  S'J  =  mm[Az(UJ),A,(Ufloe)], 

„  _  PlVnljVnl  -  Si)  ~  PEVnEjVnE  ~  SE)  +  Pi  ~  PE  +  (jBjj  -  |B||)/2 

Pi(vni  —  Si)  —  pE(vnE  —  SE) 

,  SE  =  min[Xm( UE),Am(URoe)]. 

(18) 

Here  Ai(URoe)  and  Am(URoe)  are  the  smallest  and  largest  eigenvalues  re¬ 
spectively,  of  the  Roe-averaged  matrix.  Roe-averaged  matrix  A(U),  which 
satisfies  the  following  property: 

F(Ue)  -  F(U/)  =  A(U)(Ue  -  Uj)  (19) 

Correspondingly,  A/(U/)  and  Am(U E)  are  the  smallest  and  largest  eigen¬ 
values  of  the  left  and  right  states  of  the  matrix  Ap.  The  positivity  of  pressure 
and  density  using  the  MHD-HLLC  flux  has  been  demonstrated  in  [25].  Ana¬ 
lytic  results  have  shown  that  the  flux  resolves  isolated  contact  discontinuities 
and  fast  waves  accurately. 

2.4  Implementation  of  the  Viscous  Terms 

The  viscous  terms  are  evaluated  in  two  steps.  First,  we  obtain  the  spatial 
derivatives  of  the  primitive  variables  using  the  discontinuous  Galerkin  ap¬ 
proach.  Then,  we  repeat  the  process  for  each  of  the  viscous  fluxes  using 
these  derivatives.  If  we  employ  Dirichlet  boundary  conditions  for  the  mo¬ 
mentum  and  energy  variables,  we  set  these  terms  explicitly  after  the  fluxes 
have  been  evaluated  and  then  project  the  result  using  the  orthogonal  basis. 
Here  we  use  the  average  of  the  variables  and  fluxes  at  the  interface.  This 
approach  leads  to  sub-optimal  performance  at  low  polynomial  order  p.  How¬ 
ever  it  does  not  make  much  difference  at  high  polynomial  order  p.  For  more 
details,  see  [16]. 
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3  Convergence  and  Simulations 


3.1  Convergence  Test 

A  simple  test  for  the  compressible  MHD  component  of  the  algorithm  we 
developed  is  to  consider  a  steady  irrotational  magnetic  field  and  zero  velocity. 
The  test  was  performed  as  an  initial  value  problem,  and  the  following  exact 
solution 


p  =  l,u  =  0,v  =  0 
A-iny) 

A  =  19.84+ — 2 — ,  Ee  =  9.92 

Bx  =  —cos(nx)e('~ny\  By  =  sin(nx)e^~ny^ 

was  used  as  the  boundary  conditions  and  as  the  initial  condition.  This  so¬ 
lution  but  without  the  part  concerning  Ee  was  derived  by  Priest  [29].  The 
irrotational  magnetic  field  implies  that  the  Lorentz  force  is  zero  so  the  mo¬ 
mentum  equations  are  trivially  satisfied.  The  magneto-viscous  term  is  zero 
and  the  v  x  B  term  is  also  zero.  Thus,  the  compressible  single-fluid/two- 
temperature  MHD  equations  are  satisfied.  Here,  we  take  rnt  =  1836me  [30] 
and  we  set  the  interaction  term  $  =  0  to  make  the  problem  simpler  (m,  is 
the  ion  mass  and  me  is  the  electron  mass). 

The  domain  and  hybrid  discretization  we  used  are  depicted  in  figure  2. 
We  also  show  that  the  approximation  error  decreases  exponentially  with  in¬ 
creasing  expansion  order  for  all  the  three  forms  of  error  considered  in  the  L ^ 
,  Hi,  and  +2  norms. 


3.2  Numerical  Tests  for  the  MHD-HLLC  Interface  Flux 

To  verify  our  two-dimensional  discontinuous  Galerkin  solver  with  the  MHD- 
HLLC  flux,  we  use  the  one-dimensional  benchmark  MHD  shock-tube  problem 
developed  by  Brio  and  Wu  [31].  The  one-dimensional  Riemann  problem  is 
given  for  x  G  [—1,1]: 


f  (1.000,  0,0,0, +1,0, 1.0)  for  x  <  0 

Up  (p,  Ux,  Uy ,  UZ ,  By,  Bz,  p)  \ 

[  (0.125,0,0,0,-1,0,0.1)  for  x  >  0 

(20) 
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Figure  2:  Magneto-hydrostatic  test  case.  Top-Left:  spectral  element  hybrid 

mesh  consisting  of  38  triangular  elements  and  22  rectangular  elements.  Top- 

Right:  Magnetic  streamlines  of  steady  solution;  Bottom:  Dependence  of 

steady  state  log- 

error  of  magnetic  fieltf  Bx  on  expansion  spectral  order. 

with  Bx  =  0.75  and  7  =  2.  The  solution  at  t  —  0.2  is  shown  in  figure  3,  which 
includes  the  left-moving  waves:  the  fast  rarefaction  wave,  the  intermediate 
shock  attached  by  a  slow  rarefaction  wave;  and,  the  right  moving  waves: 
the  contact  discontinuity,  a  slow  shock,  and  a  fast  rarefaction  wave.  The 
results  are  obtained  by  the  two-dimensional  discontinuous  Galerkin  solver  on 
a  mesh,  consisting  of  800  square  elements.  The  solid  line  is  the  result  using 
the  MHD-HLLC  flux  ,see  equation  (16),  the  dashed  line  is  obtained  using 
the  Lax- Friedrichs  flux,  see  equation  (15).  We  can  see  the  MHD-HLLC  flux 
gives  much  sharper  resolution,  especially  at  the  contact  interface. 

The  second  Riemann  problem  is  given  by: 

f  (1.000,  0,0,0, +1,0, 1000.0)  for  x  <  0 

Up  (/9,  "Uj;,  Up,  "U2,  Hy,  p)  \ 

[(0.125,0,0,0,-1,0,0.1)  for  a;  >  0 

(21) 

with  Bx  =  0  and  7  =  2.  This  problem  is  used  to  evaluate  the  code  for 
high  Mach  number  flow.  If  one  regards  the  term  p  +  |  \B\ 2  as  the  “hydrody¬ 
namic  pressure”,  the  system  becomes  a  standard  hydrodynamical  Riemann 
problem.  The  computational  domain  is  taken  to  be  [-1,1]  with  400  square  el¬ 
ements.  The  solution  at  t  —  0.012  is  shown  in  figure  4,  which  shows  that  the 
MHD-HLLC  flux  can  resolve  the  high  Mach  number  waves  more  accurately 
than  the  Lax-Friedrichs  flux. 

3.3  Flow  Past  a  Cylinder 

Next,  we  consider  the  problem  of  plasma  flow  with  uniform  free  stream  prop¬ 
erties  past  a  circular  cylinder.  As  the  mass  and  thermal  properties  for  elec¬ 
trons  and  ions  are  quite  different,  at  the  final  steady  state,  they  will  have 
quite  different  temperature  distributions  around  the  cylinder.  Here  we  con¬ 
sider  electrons  and  ions  having  the  same  temperature  over  the  cylinder  sur¬ 
face  as  the  free  stream  temperature.  To  simplify  our  calculation,  we  set  the 
atomic  number  of  ions  Z  —  1  and  m,  =  1836me.  We  perform  simulations  us¬ 
ing  unstructured  meshes  for  all  the  subsonic  and  transonic  cases;  it  is  shown 
in  figure  5  top,  consisting  of  490  triangular  elements;  the  mesh  used  for  the 
supersonic  case  is  shown  in  figure  5  bottom,  consisting  of  1132  triangular 
elements. 

h-refinement  is  employed  around  the  shock,  based  on  the  following  crite- 
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(a) 


(b) 


(c) 


(d) 


(e) 
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Figure  3:  ID  coplanar  MHD  Riemann  problem:  Solid  line:  MHD-HLLC  flux; 
Dashed  line:  Lax-Friedrichs  flux,  (a):  p\  (b):  p\  (c):  Up  (d):  Uy\  (e):  By. 


(a) 


(b) 


Figure  4:  High  Mach  number  Riemann  problem:  Solid  line:  MHD-HLLC 
flux;  Dashed  line:  Lax-Friedrichs  flux,  (a):  p;  (b):  p;  (c):  By ;  (d):  Ux. 
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Figure  5:  Top  :  Computational  domain  for  plasma  flow  past  a  circular  cylin¬ 
der  simulations  at  Mach  number  0.7  and  Reynolds  number  100.  Top-Left  : 
Entire  domain.  Top-Right  :  Zoom  around  the  cylinder.  Bottom  :  Computa¬ 
tional  domain  at  Mach  number  2  and  Reynolds  number  100.  Bottom-Left  : 
Entire  domain.  Bottom-Right  :  Zoom  around  the  cylinder. 
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Figure  6:  Left:  Variable  order  computational  domain;  Right:  Variable  order 
around  the  cylinder. 

rion  (the  gradient  of  density  in  the  direction  of  the  flow): 

V(p)  ■  V  >  M;  (22) 

where  M  is  an  adjustable  parameter.  The  elements  chosen  are  split  into 
4  smaller  elements.  Figure  6  plots  the  computational  domain  showing  the 
variable  polynomial  order;  p-refinement  is  used  away  from  the  shock.  To 
preserve  solution  monotonicity,  we  lower  the  spectral  order  around  discon¬ 
tinuities  appropriately.  The  polynomial  order  p  is  determined  based  on  the 
area  of  the  element  but  more  sophisticated  criteria  can  be  used  to  find  the 
optimal  p.  Due  to  the  stability  issue,  first-order  schemes  or  limiters  have  to 
be  implemented  around  the  shock.  Since  limiters  will  give  more  smearing 
results,  we  prefer  to  use  low  order  elements  in  conjunction  with  h- refinement 
close  to  the  shock. 

In  all  the  cases,  we  set  Re  =  100,  Bx  =  0.1  and  By  =  0.0  at  the  inflow, 
where  Bx  is  the  ^-component  of  the  magnetic  field  and  By  the  ^/-component 
of  the  magnetic  field.  The  two  simulations  were  run  with  polynomial  order 
p  =  5  until  the  flow  reaches  a  time-periodic  or  a  steady  state;  p-refinement 
tests  have  shown  only  very  small  difference  in  the  results. 

In  figures  7  and  8,  we  plot  the  ion  and  the  electron  temperature  contours  of 
the  instantaneous  field.  A  von  Karman  vortex  street  develops  in  the  subsonic 
regime.  However,  the  flow  is  steady  at  supersonic  states  as  shown  in  figure 
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Figure  7:  Instantaneous  nondimensional  temperature  contours  for  flow  past 
a  circular  cylinder  at  Mach  number  0.7,  Reynolds  number  100  and  cylinder 
wall  temperature  =  Te  =  1.8367.  Left  :  Ion  temperature  contours.  Right 
:  Electron  temperature  contours. 
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Figure  8:  Instantaneous  nondimensional  temperature  contours  for  flow  past 
a  circular  cylinder  at  Mach  number  2,  Reynolds  number  100  and  cylinder 
wall  temperature  Tt  =  Te  =  0.225.  Left  :  Ion  temperature  contours.  Right  : 
Electron  temperature  contours. 
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8.  In  particular,  we  have  performed  several  simulations  for  the  compressible 
flow,  using  the  one-temperature  model  and  the  two-temperature  model  in 
order  to  present  differences  in  the  forces  and  frequency. 

Specifically,  in  order  to  compare  the  one-temperature  model  with  the 
two-temperature  model,  we  integrate  the  flow  field  over  one  time  period 
and  obtain  the  time- averaged  flow  held.  In  figures  9  and  10,  we  plot  the 
time-averaged  contour  lines  of  the  ion  and  the  electron  temperature  using 
the  two-temperature  model  and  the  temperature  contours  using  the  one- 
temperature  model.  Also,  we  present  temperature  profiles  along  a  line  on 
top  of  the  cylinder  aligned  with  the  vertical  axis,  as  shown  in  figure  11,  to 
compare  the  temperature  distributions.  In  figures  12  and  13,  we  plot  the 
normal  profiles  of  temperature  in  one-temperature  model  and  profiles  of  the 
electron  and  the  ion  temperature  in  two-temperature  model  starting  at  the 
top  of  the  cylinder  at  Mach  number  0.7  and  2.  From  these  figures,  we  can 
see  that  the  results  of  one-temperature  model  are  quite  different  from  the 
ones  of  the  two-temperature  model.  In  figures  12  and  13,  electrons  and 
ions  have  the  same  temperature  on  the  cylinder  surface  and  the  electron 
temperature  increases  while  the  ion  temperature  decreases  in  the  direction 
away  from  the  cylinder  surface.  From  figures  9  and  12,  we  can  see  that 
the  ion  temperature  is  larger  than  the  electron  temperature  ahead  of  the 
cylinder  at  Mach  number  0.7.  However  the  electron  temperature  is  larger 
than  the  ion  temperature  on  the  two  sides  of  cylinder.  Figure  13  shows  the 
electron  temperature  is  smaller  than  the  ion  temperature  on  the  two  sides 
of  cylinder  at  Mach  number  2.  Generally, temperature  profiles  obtained  from 
one  temperature  model  is  similar  to  the  ion  temperature  rather  than  the 
electron  temperature.  The  relation  between  two  temperatures  model  and 
one  temperature  model  can  be  given  as:  T  =  TlXTe ,  if  ions  and  electrons 
have  the  same  number  density  nl  =  ne. 

4  SUMMARY 

We  have  developed  a  discontinuous  Galerkin  solver  to  model  two-temperature 
plasmas  as  part  of  a  hierarchical  modelling  approach  and  a  compromise  be¬ 
tween  single-temperature  and  two-fluids  models.  We  have  demonstrated 
spectral  convergence  for  an  analytical  problem  and  also  demonstrated  the 
robustness  of  the  method  in  dealing  with  shocks  without  the  use  of  flux 
limiters  or  artificial  viscosity  terms.  The  issue  of  preserving  positivity  is 
addressed  by  introducing  a  new  interface  flux,  the  MHD-HLLC  flux. 


22 


Figure  9:  Left  :  Time  average  contour  lines  of  ion  nondimensional  temper¬ 
ature  (upper  half  plane)  and  electron  nondimensional  temperature  (lower 
half  plane)  at  Mach  number  0.7,  Reynolds  number  100  and  cylinder  wall 
temperature  T)  =  Te  =  1.8367  using  the  two-temperature  model.  Right  : 
Time  average  contour  lines  of  nondimensional  temperature  using  the  one- 
temperature  model. 
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Figure  10:  Left  :  Time  average  contour  lines  of  ion  nondimensional  tempera¬ 
ture  (upper  half  plane)  and  electron  nondimensional  temperature  (lower  half 
plane)  at  Mach  number  2,  Reynolds  number  100  and  cylinder  wall  tempera¬ 
ture  Ti  =  Te  =  0.225  using  the  two-temperature  model;  Right  :  Time  aver¬ 
age  contour  lines  of  nondimensional  temperature  using  the  one-temperature 
model. 


Figure  11:  Location  where  profiles  of  temperature  along  the  line  shown 
aligned  with  the  vertical  axis  are  taken.  Contours  of  ion  temperature  at 
Mach  number  0.7  are  shown  in  the  background. 
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Figure  12:  Normal  profiles  of  temperature  nondimensionalized  by  cylinder 
wall  temperature  at  Mach  number  0.7,  Reynolds  number  100  and  cylin¬ 
der  wall  temperature  T%  —  Te  —  1.8367.  Solid  line:  Ion  temperature 
from  two-temperature  model;  Dotted  line:  Electron  temperature  from  two- 
temperature  model;  Dashed  line:  Temperature  from  one-temperature  model. 


Figure  13:  Normal  profiles  of  temperature  nondimensionalized  by  cylin¬ 
der  wall  temperature  at  Mach  number  2,  Reynolds  number  100  and  cylin¬ 
der  wall  temperature  T*  =  Te  =  0.225.  Solid  line:  Ion  temperature 
from  two-temperature  model;  Dotted  line:  Electron  temperature  from  two- 
temperature  model;  Dashed  line:  Temperature  from  one-temperature  model. 
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Examination  of  contours  of  the  divergence  of  the  magnetic  held  in  how 
past  a  cylinder  revealed  some  non-zero  values  in  the  region  around  the  rear 
stagnation  point.  However,  this  was  not  growing  in  time  and  thus  no  numer¬ 
ical  instabilities  were  induced  even  after  long-time  integration.  This  diver¬ 
gence  held  can  be  totally  eliminated  by  occasional  Helmholtz  decomposition 
of  the  magnetic  hux  vector,  see  [16]  and  [26].  We  are  currently  developing 
a  discontinuous  Galerkin  method  for  two-huids  plasmas  and  we  will  report 
those  results  in  a  future  publication. 

APPENDICES 

A:  Detailed  flux  terms  in  single-fluid/two-temperature 
equations 

We  present  here  in  detail  the  hux  terms  involved  in  the  governing  equations 
of  the  single-huid/two-temperature  plasma.  Many  of  the  parameters  used 
are  listed  in  table  1. 


^ I deal 


^ I deal 


y 


? I deal 


i Vise 


■\Visc 


( pu ,  pu2  —  B2  +  p,  puv  —  BxBy,  puw  —  BXBZ,  0,  uBy  —  vBx,  uBz  —  wBx, 

( E  +  p)u  -  (v  •  B )BX,  (ee  +  pe)u)T 

(pv,  pvu  —  ByBx,  pv2  —  B2  +  p,  pvw  —  ByBz,vBx  —  uBy,  0,  vBz  —  wBy , 

(E  +  p)v  -  (v  •  B )By,  (ec  +  Pe)v)T 

(pw,  pwu  —  BZBX,  pwv  —  BzBy,  pw2  —  B2  +  p,  wBx  —  uBz,  wBy  —  vBz, 

0,  (E  +  p)w  -  ( v  •  B )BZ,  (ee  +  pe)w)T 

2  .  du  1  1  .du  dv ,  1  .  du  dw ,  1  .  dBv  d  Bx. 

wfv  -  W -). yW  +  &>•  ydy  +  &  >■  °- -y 1 -  ar)' 


Sm  dx 
1  ,dB, 


3 

f)B 


Sr '  dx 

1  1  <9(|B|2 


1  ,  2 . _  .  „  l<9v2  1  dE  1  dTe 

dTh  yd_3(v  'v>“  +  v '  v“  +  2  +  W, H  +  pF.  dV 


Sr  v  2  dx 


B  -VBX), 


dT, 


e \T 


SvePr e  dx 


dB„ 


1  dv  du  2  dv  1  1  dv  dw  1  ,dBx 

+  dy>'^dy  ~~  3  V  '  +  ~dy^~~S,S~dy  ~  ~dx 


),0, 


1  ,dBr  dB, 


V  dy 


1.2  .  „  1  <9v2  1  dE  1  dTe 

dz  SVi  3  ^  ;  2  dy  Pr i  dy  Pre  dy 
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Variable 

Description 

p(x,  t)  =  pi  +  pe 

single  fluid  density 

v(x,t)  =  (u,v,w)(x.,t) 

single  fluid  velocity 

B(x,  t)  =  (Bx,By,Bz)(x,t) 

magnetic  fields 

E  =  Bv  +  §(pv  '  v  +  B  '  B) 

total  energy 

P  =  Pi  +  Pe 

total  pressure 

P  =  P  +  z B  •  B 

pressure  plus  magnetic  pressure 

_  Pi  y  _  Pe 

1  Pi'1  6  Pe 

ion  and  electron  non-dimensional  Temperature 

Pr.  —  cpPi  pr  —  CpPe 

1  '  i  K.  i  1  '  e  K 

rv7.  rve 

ion  and  electron  Prandtl  Number 

V 

Magnetic  resistivity 

l^ii  l^e 

ion  and  electron  Viscosity 

c  _  poVa^o  c  _  poVa^o 

m  Pi  ’  ve  Me 

ion  and  electron  Viscous  Lundquist  number 

C  _  P-oVaX 0 

Dr  —  „ 

Resistive  Lundquist  number 

Cp 

Specific  heat  at  constant  pressure 

vj  =  ^ 

A  Pop 

Alfven  wave  speed 

A  - 

Alfven  Number 

Table  1:  Variables  and  parameters  used  in  the  equations  of  single-fluid/two- 
temperature  compressible  MHD. 


i Vise 


s 


MHD 


1  ,10(|B| 


Sr  2  dy 

1  .dw 

du 

’  Svi  ^  dx  dz 

1  (dBy 

dBz 

SB  dz 

dy 

1  13(|B| 

2)_ 

—  B  ■  VB. 
1 


d  Tf 


e\T 


vh  SvePre  dy  ' 
dw  dv.  2  ,  dw 


dy  dz  ’  Sm  dz 


1, 

3 


1  ,dB,r  dB, 


SB  dz 


dx 


1  dw2 


)'0'S-(“3(V’V)”’  +  V’V“’+2  'dz 


+ 


j_dT] 

Pr,  dz 


+ 


\_dj\ 

Pre  dz 


—  B  ■  VB 


dT( 


e \T 


SvePre  dz 

1  _  1 


=  (0,  0,  0,  0,  0, 0,  0,  0,  (v  •  X7pe  +  —  re  :  Vv  +  —  (V  x  B)  ■  (V  x  B ))) 


5, 


Sr 
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B:  From  two-fluid  equations  to  one-fluid  two-temperature 
equations 


The  magnetohydrodynamic  model  treats  the  plasma  as  a  single  fluid.  In  the 
following,  we  derive  the  single- fluid/two-temperature  plasma  equations  from 
the  two  fluids  plasma  equations.  We  denote  the  electron  and  ion  masses  by 
m  and  M  respectively,  rq  and  ne  are  the  ion  and  electron  number  density. 
ji,  je,  Ti ,  Te,  pi,  pe  and  Ei,  Ee  are  the  current  density,  temperature  pressure 
and  hydrodynamic  energy  density  of  ion  and  electron  respectively.  Etot  and 
ee  are  the  total  energy  density  and  electron  internal  energy  density. 

We  Define: 


M  = 

P  = 
ji  = 

je 

,i  = 
P  = 
T  = 

v  = 


Ei 

Ee 

Etot 


1836rq  rq  =  ne  =  n 

riiM  +  nem  =  n(M  +  m)  ~  nM 

newi 

—neve 

ne(vi  -  ve)  =  ne(j*  +  je) 

Pi  +  Pe 

Ti  +  Te 


n(M  v*  +  mvc 


Mwi 


P 

-  mv. 


M  +  m 

mj 


m 

V‘  +  MV' 


Mne 

j  j 

y  _  _  _ 

ne  ne 

nkTi  1  2 

=  ^i  +  2nMv‘ 

nkTe  1  2 

=  — +  -mnve 

—  Ei  +  Es  +  — B' 


He 


i  /  . 

7  -  1 


where  7  —  5/3  and  k  is  Boltzmann  constant. 


(23) 

(24) 

(25) 

(26) 

(27) 

(28) 

(29) 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 

(36) 


We  make  following  assumptions: 

1.  Quasineutral  approximation. 

n  =  7ij  =  ne  (37) 


2.  High  collisionality.  We  assume  both  the  electrons  and  ions  are  collision 
dominated.  The  collisions  rapidly  randomize  the  distribution  function  giving 
rise  to  an  isotropic  pressure. 

3.  Small  Larmor  radius. 

Generalized  Ohm’s  Law: 


E 


V j  +  v  x  B 


j  X  B  -  Vpe 


ne 


(38) 


With  the  ‘small  Larmor  radius’  approximation,  the  generalized  Ohm’s 
law  is  simplified  as: 

E  =  ?|j-vxB  (39) 

4.  Electrons  move  much  faster  than  ions. 

| vi  |  <  |ve|  (40) 


Mass  conservation 

In  two  fluids  plasma  mass  conservation  equations,  we  have: 

0  (41) 

0  (42) 

By  Multiplying  the  ion  and  electron  masses  M  and  m,  respectively,  and 
adding  above  two  equations  together,  we  produce  the  ‘single-fluid  mass  con¬ 
servation  equation’, 

dU^Mdt  ^  +  V  ’  +  mVe)l  =  fgf  +  V  ’  (dv)  =  0  (43) 


UUi  ,n  t  \ 

—  +  V  •  (mvi)  = 

dne  „  , 

—  +  V  •  (neve)  = 
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Momentum  conservation 

In  two-fluids  plasma,  the  momentum  conservation  equations  are: 


d(nMv,)  +  v  (nMv.vty  =  _ Vp.+ne( E  +  v,xB) 
+V  •  (fjliTi) 

d  n _|_  y  .  (nmveVg!)  =  —  Vpe  —  ne( E  +  ve  x  B) 


(44) 


dt 

+V  •  (/xerc)  (45) 

ere  r*  =  (djVU  +  dpvij)  -  |V  •  and  re  =  (djVei  +  <9*vei)  -  |V  •  ve4?;j. 
By  adding  the  above  two  equations  together,  we  have: 

+  V  •  (nMvjV*I  +  nmveVgI)  «  +  V  •  (pvv*I) 

=  -V(p*  +  pe)  +  ne(vj  -  ve)  x  B  +  V  ■  (/q7*  +  peTe) 

=  -  Vp  +  j  X  B  +  V  •  (HiTi  +  peTe) 


(46) 


where  we  define  p  =  p*  +  pe  and  j  =  ne(v,  —  ve). 

Since  j  x  B  =  —(V  x  B)  x  B  =  ——V  ■  (— BB*  +  ^ IB |2I) ,  the  combined 
’single-fluid  momentum  conservation  equation’  is  obtained  as  following, 

=  -V  •  (pW  -  +  (p  +  -^|B|2)I  -  Pin  -  peTe)  (47) 

at  un  2pa 


After  non- 


dimensionalization,  we  have: 

-  BB'  +  (p+  ^IBI2)!  -  J-Ti  -  ^-Te) 


<9(pv) 

-dT  =  ~v'  (',vv 


■'ve 


where  SA  =  and 

fJ'i  l^e 


(48) 


Magnetic  Field 


<9B 

lit7 


=  -V  x  E 


where  E  =  —V7  X  B  —  v  x  B. 

0  Mo 


<9B 


=  —V  x  (B  x  v 


-V  x  B) 


(49) 
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(50) 


where  v  =  iUY}tmVe  •  After  non-dimensionalization,  we  have, 


<9B  1 

—  = -V  x  (B  x  v+ — V  x  B)  (51) 

where  Sr  =  PoVaLo/t]. 

Energy  conservation 

In  two-fluids  plasma  energy  conservation  equations,  we  have: 


dEt 

dt 

dEe 

dt 


+  V  •  [(Ei  +  Pi)vi\  =  jj  •  E  +  V  •  ■  Ti  +  KiVTi)  (52) 

+  V  •  [{Ee  +  pe)ve]  =  je  ■  E  +  V  ■  (/ieVe  •  Te  +  Ke VTe)  (53) 


By  substituting  V;  m  v  +  and  ve  «  v  —  -J-  into  the  above  two  equations, 
we  have: 


9Ei  +  V  •  [(E,;  +  Pi)v]  +  ■  [(Ei  +  Pi)(V  x  B)]  = 


at 

<9E, 


M  /jL0ne 
•  E  +  V  •  (niVi  ■  Y  +  KiVTi ) 


+  V  •  [(Ee  +  pe)v] - — V  •  [(Ee  +  pe) (V  x  B)]  = 

C'L'  I*Jj  qTI/G' 

je  ■  E  +  V  ■  (fteVe  ■  Te  +  KeVTe) 


(54) 


(55) 


Adding  the  above  two  equations  and  considering  in  two-dimensional  space,  V- 
(V  x  B)  =  0,  the  above  equation  can  be  simplified  as: 


d  E 


tot 


at 


+  V  ■  [(Etot  +  p)v]  =  j  E 
+  V  ■  (p,iVi  ■  Ti  +  pewe  ■  Te  +  ki VTi  +  fi)eVTe) 


(56) 


Using  the  magnetic  held  equation  obtained  above,  we  have: 
1  <9B^ 

=  — B  •  V  x  (B  x  v  +  r/V  x  B) 


2  dt 


-  —  V  •  [(-|B|2/  —  BB*)  •  v  +  r/j  x  B]  —  j  ■  E 


(57) 
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Adding  the  above  two  equations,  we  have: 


-QjT-  +  V  •  [(Etot  +  p)v]  =  j  •  E  +  V  •  (fjliVi  ■  Ti  +  [iewe  ■  Te  +  kiWTi  +  ke\7Te ) 

-  V  ■  [(^|B|2/  —  BBf)  •  v  +  r/j  x  B]  —  j  •  E  (58) 

Since 

|B|2 

(V  x  B)  x  B  =  B  •  VB  -  V^-  (59) 

We  have  the  combined  energy  conservation  equation: 

dE  1 

—EL  =  -V  •  [(Etot  +p)v  +  (-|B|2/  -  BB')  •  v  -  fjtiv  •  t£  -  /xev  •  re 

-  fcivri-fcevre  +  77( b-vb-vK)]  (eo) 

After  non-dimensionalization,  we  have, 


-v  ■  [(Em+p)v  +  (5b|V  -  BB1)  •  v  -  — j— VT,  -  Vre 

Z/  ovi±  i  2  •  e 

^v-re  +  ^(B-  VB-  V®-)]  (61) 


where  Pry  =  and  Pre  =  — . 

Avj  Ave 


Electron  energy  conservation 

The  electron  energy  conservation  equation  is: 
dc 

-q y  +  V  •  [(ee  +  pe)ve]  =  V  •  (ne VTe)  +  /iere  :  Vve  +  ve  ■  Vpe  +  r?j  •  j  (62) 


By  substituting  ve  k  v  -  ^  into  above  equation,  we  have: 

f)f  1 

—A  +  V  ■  [(ee  +  pe)v] - —  V  •  [(ee +  pe)(V  x  B)] 

Cji/  /jl0tlc 

=  V  •  (KeVPe)  +  /xere  :  Vve  +  v  ■  Vpe  -  X  •  Vpe  +  r/j  •  |63) 

p0ne 
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Considering  the  two-dimensional  space.  V  •  (V  x  B)  =  0  and  (-VxB'>  ■  Vne  =  0, 
above  equation  can  be  simplified: 

^  +  V  •  [(ee  +  pe)v] 

=  V  •  (Kcvre)  +  /deTe  :  Vv  +  V  •  Vpe  +  rj]  •  j  (64) 

After  non-dimensionalization,  we  have, 


-V  •  [(ee  +  pe)v  -  1  vre] 

^ve-L  '  e 

v  ■  Vpe  +  -^-re  :  Vv  +  -^-(V  X  B)  •  (V  X  B)  (65) 
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